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We perform Monte Carlo simulations of large two-dimensional Gaussian Ising spin glasses down to 
very low temperatures (5 = 1/T = 50. Equilibration is ensured by using a cluster algorithm including 
Monte Carlo moves consisting of flipping fundamental excitations. We study the thermodynamic 
behavior using the Binder cumulant, the spin-glass susceptibility, the distribution of overlaps, the 
overlap with the ground state and the specific heat. We confirm that T c — 0. All results are 
compatible with an algebraic divergence of the correlation length with an exponent v. We find 
— X/v — —0.295(30), which is compatible with the value for the domain-wall and droplet exponent 
6 ~ —0.29 found previously in ground-state studies. Hence the thermodynamic behavior of this 
model seems to be governed by one single exponent. 

PACS numbers: 75.50.Lk, 05.70.Jk, 75.40.Mg, 77.80.Bh 



I. INTRODUCTION 

Spin glasses^ are the prototype models for disordered 
systems investigated extensively in statistical physics 
during the last three decades. These systems exhibit 
complex energy landscapes resulting in many interest- 
ing phenomena, like glassy behavior and aging. Despite 
much effort, many questions are still open. 

For two-dimensional spin glasses with only nearest- 
neighbor interactions, it is now clear that no stable spin- 
glass phase at finite temperature existsS* 3 ^*^. This 
means that they are paramagnetic at any finite tem- 
perature and that they exhibit spin-glass behavior only 
at T = 0. Thus, it is widely believed, that their 
low-temperature behavior can well be described by the 
droplet theoryiSiS. The droplet picture assumes that the 
low-temperature behavior is governed by droplet-like ex- 
citations, where excitations of linear spatial extent I typ- 
ically cost an energy of order l e . Thus in the thermody- 
namic limit the excitations which flip a finite fraction of 
the spins cost an infinite amount of energy if 9 > 0. 
For the two-dimensional spin glass with Gaussian in- 
teractions, since it exhibits no order at T > 0, 9 < 
holds— Furthermore it is usually assumed that the en- 
ergy of different types of excitations, e.g. droplets and 
domain walls, induced by changing the boundary condi- 
tions, are described by the same exponent 9. Indeed, re- 
cently is has been confirmed by calculating exact ground 
statesii that 9 = —0.282(2) for droplet and domain-wall 
excitation a 12 : 13 . For domain walls, small sizes are suf- 
ficient to see the asymptotics, hence similar values have 
been found previouslySii^ii^iSiiiiSiiS in this case. On the 
other hand, for some droplet excitations, this behavior is 
only visible for not-too small system sizes L > 50, which 
explains why in a similar preceding studySS of smaller 
sizes an apparently different exponent close to 9 = —0.47 
has been found. For other types of droplet-like excita- 
tions, the exponent 9 rj —0.28 is again already visible for 



small sizes 21,22 . Hence, the behavior at zero temperature 
seems to be relatively clear—. 

The situation is different for the small but finite- 
temperature behavior. In case the correlation length £ 
diverges algebraically for T — > like £ ~ T~ v , the criti- 
cal exponent v of the correlation length is related through 
a simple renormalization argumenljiaii^ to the droplet ex- 
ponent 9 via 9 = —1/v. Several numerical studies to ob- 
tain v at finite temperatures have been performed. Using 
transfer-matrix calculations of long (L x up to 10 6 ) and 
narrow (L y up to 11) stripes, values ofi& v — 2.96(22) 
respectively 2 ? v — 4.2(5) have been found. For small 
(L = 10) square systems^, a value of v = 2.1(1) was 
found. Also a couple of Monte Carlo (MC) simulations 
have been performed. For small sizes (L = 12) and rela- 
tive large temperatures T > 1, a value of v = 3.6(2) has 
been founds. Later, for similar system sizes but lower 
temperatures, a value of v = 1.8(4) was obtained— Fur- 
thermore, a cluster Monte Carlo simulation of three large 
samples (L — 128) in the temperature range T G [0.4, 1] 
has bee performed, resulting 3 ° in v — 2.0(2). 

Since the finite-temperature studies performed so far 
suffer from either too small systems or too large tem- 
peratures (or both), we have performed a Monte-Carlo 
study of large systems up to size L = 75. By applying 
a recently developed cluster algorithm^, in connection 
with an extension presented below, we are able to equi- 
librate the system at much smaller temperatures than it 
was possible before. Parallel and independently of us, 
H.G. Katzgraber, L.W. Lee and A. P. Young performed 
a related study 3 ^ using a similar algorithm. They focus 
on the direct calculation of the correlation length, while 
we study here other thermodynamic quantities like the 
Binder parameter, the spin-glass susceptibility, the dis- 
tribution of overlaps, the overlap with the ground state 
and the specific heat, and we infer from these results the 
asymptotic behavior of the correlation length indirectly 

The model we study consists of N = L 2 Ising spins 
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Si = ±1 on a square lattice with the Hamiltonian 

H=-^2J ii S i S j , (1) 

(id) 

where the sum runs over all pairs of nearest neighbors 
and the Jij are quenched random variables which 
are distributed according to a Gaussian distribution with 
zero mean and unit variance. In the following, we denote 
the thermal average by (...) and the quenched-disorder 
average by [. . Periodic boundary conditions in both 
directions are applied. 

The rest of the paper is organized as follows. First we 
explain the algorithm we have applied and show that it 
is able to equilibrate the system. In the main part, we 
present our results for the different thermodynamic quan- 
tities mentioned above. In the last section, we summarize 
our results. 



II. ALGORITHM 

We have made extensive MC simulations of our system. 
To reach equilibrium down to very low temperatures for 
large system sizes, we have used a recently developed 
cluster algorithm (details can be found elsewhere^). To 
speed up equilibration at very low temperatures we have 
devised and used a new procedure. It consists in main- 
taining a list of the lowest-energy elementary excitations, 
and to use flipping them as MC moves. This procedure 
works as follows: 

1. Compute the ground state of the system. To do this 
we used a heuristic renormalization-group based 
algorithmic^ To test the efficiency of the method, 
we checked, for systems of size L = 100 with open 
boundaries, that this algorithm systematically finds 
the true ground state by comparing with the result 
of an exact algorithm^ (which works only for pla- 
nar graphs, i.e. not with fully periodic boundary 
conditions). Moreover during the production runs 
using the MC simulations, we never found excita- 
tions with negative energy, which confirms that the 
true ground state has been found. 

2. During the equilibration phase, we systematically 
compare the low temperature (T < 0.2) configura- 
tions with the ground state. They differ by clusters 
that are flipped. We maintain a list of the lowest- 
energy excitations thus found (we keep up to 10000 
of such excitations for the largest system). Note 
that we consider only elementary excitations whose 
boundary is connected (two clusters flipped inside 
one another define two independent elementary ex- 
citations, not one). Then we introduce (in addition 
to the cluster algorithm) a new Monte Carlo move 
at low temperature (T < 0.2): choose an elemen- 
tary excitation in the list and try to flip it using 
the Metropolis criterion. As soon as the list no 



-6000 



E 




MC steps 

FIG. 1: Sample equilibration test for L = 75 and /3 = 5 
averaged over 200 samples, the curves (Eq.|5J start to overlap 
around the equilibration time (here 25000 MC steps). 



longer evolves (i.e. when all the first excitations 
have been found) these moves trivially respect de- 
tailed balance. 

3. During the production phase, we no longer try to 
find new excitations, we simply use the cluster al- 
gorithm together with single-spin flips. 

Note that we obtain a list of the first elementary exci- 
tations as a by-product of this algorithm. To check that 
equilibrium has been reached, we have used the criterion 
described in R.ef. lsil which is based on the following iden- 
tity valid for a Gaussian distribution of the interactions: 



[E]j = - (2) 

where E = (hi) is the average energy, = 1/T the in- 
verse temperature, TV; = 2N is the number of bonds, and 
qi denotes the link overlap between two independently 
chosen configurations {S"}, {Sf} for the same disorder: 



(id) 

In Fig.Q]we show how both sides of Eq. 0evolve during 
equilibration. The system can be considered equilibrated 
when the curves start to overlap. Note that within our 
algorithm, the temperature, where equilibration takes in 
longest time, is not the lowest temperature, because the 
excitation flips are most efficient at the lowest tempera- 
tures. For example at the lowest temperature, the ground 
state may be found after the very first step of the al- 
gorithm because all excitations present in the starting 
configurations can be flipped at once. 
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FIG. 2: Binder cumulant g as a function of inverse tempera- 
ture P for different system sizes L. 



III. RESULTS 

Wc considered the following sizes for our simulations: 
L = 10, 25, 35, 50 and 75 over a very large range of 
temperatures 0.02 < T < 5. We respectively treated 
1000, 1000, 1000, 500, 200 samples for the different sizes. 
We simulated at different values of the temperature, the 
number of different values ranging from 19 to 59 (with 
increasing size). For each sample, we simulated 64 inde- 
pendent configurations at each temperature. 

To study our system, we have measured different quan- 
tities and averaged over the different samples. Our first 
quantity of interest is the Binder cumulant 36,37 g defined 
by 
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(4) 



Here q is the overlap between two independent equili- 
brated configurations {S"} and {Sf } of the same disor- 
der realization 
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(5) 



In the thermodynamic limit, the Binder cumulant is zero 
in the paramagnetic phase and around one in the spin- 
glass phase. To study g at finite sizes, we consider the 
divergence of the correlation length £ when approaching 
the transition temperature T c 



£ ~ (T — r c y 



(6) 



Since g does not show any critical behavior near the crit- 
ical point, and according to the basic assumption that 
it is a function of the relation of system size L to the 
correlation length £, g scales as: 



FIG. 3: Scaling of the Binder cumulant in the critical region 
obtained by plotting g as a function of the rescaled inverse 
temperature /3L -1 /" using 1/u — 0.295 



g~g{L/0 = g{L(T-T c y). 



(7) 



Hence, g is independent of L at T = T c which allows 
the location of T c : the curves for different system sizes L 
should intersect at T c . In Fig.[2]we show the value of g as 
a function of (3 for different system sizes L. The curves for 
g converge to a crossing point at very low temperature 
which indicate that we most probably have T c = 0, in 
accordance to the recent believeS* 3 ^^. 

As stressed in the introduction, the value of v has been 
an open question for a long time. The reason for this 
is the presence of large finite-size corrections to scaling. 
We see this, when trying to perform a finite-size scaling 
plot, i.e. when plotting g as a function of (5L~ X I V for a 
suitably chosen value of v. In fact, no value for v allows 
the whole curves to collapse on a master curve. One has 
to select only a domain of parameters near the critical 
regime, namely large L and small T (and thus large g). 
Keeping only L > 35 and g > 0.5 we find that 1/u ~ 
0.295 (v ~ 3.39). The scaling plot for this value is shown 
on Fig. |3 (note that we show all the data points, but only 
those in the range mentioned above where used to find 
the exponent): for large sizes and low temperatures, a 
very good data collapse is obtained. 

To confirm this value of v and to estimate the error, 
we need to somehow quantify the quality of the collapse 
of the curves. To do this, we use a procedure similar to 
one proposed by Kawashima and Itc 38 which we detail in 
the appendix. Using this procedure, we define a function 
S(v), measuring the quality of the fit as a function of 
the chosen value of v. S(y) should be around one if 
the collapse is good (taking the error bars into account) 
and much larger otherwise, it behaves somehow like a \ 2 
test. In Fig. 0] we show the value of S*(^) for the scaling 
according to Eq. again we used only the data close to 
T = 0. We see that the minimum corresponds to S ~ 
1.28 which is good and we can also have an estimation 
for the error on v. 0.295 ± 0.03. 

In Fig. |3 we see that the high temperature part, which 
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FIG. 5: Scaling of the Binder cumulant out of the critical 
region, when choosing 1/v = 0.38. 

was not used to obtain is, does not scale well. It is as if 
it requires another value of the exponent. In fact, if we 
use only data for which L > 25 and 0.1 < g < 0.5 we find 
1/v ~ 0.38 (and S = 6.6 which is not good) the resulting 
plot is shown in Fig. [5j The scaling is not good but it 
explains why such a high value of 1 jv appears in previous 
papers: large system sizes, low temperatures and a good 
criterion for equilibration are necessary to obtain valid 
results. 

We now turn to the spin-glass susceptibility 

X sg = L d [(q 2 )}j (8) 

(here d = 2). Since the ground state of Ising spin glasses 
with a Gaussian distribution of the interactions is unique 
(i.e. q — 1), the susceptibility shows 28 - 29 the following 
finite-size behavior at low temperatures: 

Xsg ~ L 2 x(L/0 ~ L 2 x(LTn, (9) 



FIG. 7: Scaling of the spin-glass susceptibility ;^sg 

(with i(L/£) going to a constant for T — > 0). A little 
bit away from the critical region, where the correlation 
length is small compared to the system size, the sus- 
ceptibility should not show any system-size dependence, 
hence the L 2 factor must cancel out (x(^) ~ x ~ 2 ), an d 
we obtain 

Xs S - T' 2 \ (10) 

at large L and finite T. We show our result for Xsg in Fig. 
The line added corresponds to a power law with expo- 
nent 5.0 which corresponds to l/v = 0.4. The line does 
not fit the data very well and there is some upward cur- 
vature of the data which indicates a larger value for the 
true exponent. This is fully compatible with the previous 
results, because 1/v ~ 0.295 results in 2v ~ 6.8. To test 
the prediction in Eq. we show in Fig. 0a plot of x/L 2 
as a function of flL^ 1 /" . We get a quality S = 1.64 of 
the finite-size scaling when considering data with L > 35 
and Xsg/L 2 > 0.3 which means that the scaling is good. 
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FIG. 8: The distribution P(g) of overlaps for different system 
sizes at inverse temperature (3 = 2. 



FIG. 9: The overlap g ga with the ground state as a function 
of the inverse temperature for different system sizes L. 



We also studied the full distribution P{q) of overlaps. 
The data at /3 = 2 is shown in Fig. Q3J We observe that 
the distributions are perfectly symmetric with respect to 
q = 0, which is another indication that our simulation 
is well equilibrated. Furthermore, we can see a crossover 
from a two-peak structure for small system sizes to a 
trivial distribution as the size is increased. This again is 
fully compatible with the notion that there is no ordered 
phase a low temperatures except at T = 0, hence the 
system is paramagnetic at finite temperatures. 

We furthermore looked at the overlap of the 
equilibrated configurations with the ground-state 
configuration-, this can be seen as a "spin-glass 
magnetization" 



<?gs 



(11) 



J j 



where Sf denotes one of the two ground states. This is 
expected^ to scale as 



<7 gs ~ q&(L/g) ~ q~ gs (LT u ). 



(12) 



The raw data is shown in Fig. |U1 while the data rescaled 
according to Eq. PEZlwith 1/V = 0.295 is shown in Fig.lTUl 
When considering data with L > 35 and q gs > 0.5 we find 
quality S = 2.50. The quality of this scaling is lower than 
the one for the Binder cumulant, nevertheless, the result 
still supports the findings from above. Interestingly, in 
Ref. Q a similar value 1/V = 0.28 was already found 
whereas only small systems L < 12 at high temperatures 
(f3 < 1.4) where studied. 
We also studied the specific heat 
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dT 



(13) 



Although the specific heat is not expected to show any 
singularity for the spin-glass transition, it is nevertheless 




FIG. 10: Scaling plot of the overlap with the ground state: 
as a function of /3L _1/ ", with 1/V = 0.295. 



interesting to compare with previous studies. The data 
is shown in Fig. ^2 The specific heat presents a finite 
peak around (3 = 0.82 known as the "Schottky anomaly" . 
The decay of the specific heat at small T goes roughly 
as c ~ T similarly to previous studies using transfer- 
matrix calculations on long stripes-- respectively hierar- 
chical lattices^. Nevertheless the data is not perfectly 
described by a power low even at such low temperatures 
and it is not clear what the asymptotic behavior truely 
is. 

Finally, we have also directly evaluated the correlation 
length £ by measuring spin-spin correlations as a func- 
tion of the separation of the spins and fitting suitable 
functions to the data. Since the behavior of the correla- 
tion is already studied in detail in Ref. 0, we do not go 
into details here. We only mention that our results are 
fully compatible with the results of Ref. when going 
to large system sizes and studying low temperatures, we 
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APPENDIX A: DETERMINING THE QUALITY 
OF THE SCALING LAWS 



FIG. 11: The specific heat c as a function of inverse temper- 
ature for different system sizes L. The straight line is a 
putative asymptotic behavior c ~ T when T — > 0. 



observe a behavior compatible with —1/u = —0.295. 



IV. SUMMARY AND CONCLUSIONS 

We have studied the low-temperature behavior of two- 
dimensional spin glasses with a Gaussian distribution of 
the interactions. Using a sophisticated cluster algorithm 
in combination with exploring lowest-energy excitations 
close to T — 0, we were able to equilibrate large system 
sizes up to I = 75 down to very small temperatures 
T = 0.02. 

We have studied several thermodynamical quantities, 
like Binder cumulant, susceptibility, distribution of over- 
laps, overlap with the ground state and specific heat. Our 
main findings are as follows. From the Binder cumulant, 
and the distribution of overlaps, we see that no stable 
spin-glass phase at finite temperature exists, i.e. T c = 
in accordance to recent studies. From the scaling be- 
havior of the Binder cumulant and the susceptibility, we 
find that the correlation length diverges algebraically for 
T — > 0, in contrast to the model with a bimodal distribu- 
tion of the interactions 5 . The main open question was, 
whether the exponent 9 w —0.29, describing the scaling 
of droplet and domain- wall excitations, is equal to —1/v, 
v being the exponent of the correlation length. Our re- 
sults — l/v= —0.295(30) obtained for large sizes and at 
low temperatures indeed support 9 = —\jv. 

To conclude, when taking ground-state studies into 
account: the thermodynamic behavior of the two- 
dimensional Gaussian Ising spin glass is trivial for finite 
but low temperatures and governed by one single expo- 
nent —\jv = 9 ks —0.29, as predicted by the droplet 
picture simple renormalization arguments. 



To determine the quality of the scaling laws we need 
some quality criterium that somehow measures the dis- 
tance of the data to the master curve. The difficulty 
arise from the fact that the master curve is unknown 
and must be determined from the data. Some years ago, 
Kawashima and Itc 3S proposed a method. We present 
here a refinement which, according to our experience, 
seemed to be stabler and more precise. 

After applying the scaling law, we have k sets of points. 
Each set is composed of n, points (i — 1 . . . k) of the 
form (xij , yij , dyij ) with dy being the standard error on 
y and j = 1 . . . m . In the following we suppose that 
xn < Xi2 < . . . < Xim ■ We define the quality as 
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Y 



M^rf dy 



(Al) 



where Yy and dYij are the estimated position and stan- 
dard error of the master curve at Xij . M is the number of 
terms in the sum (we only consider the terms for which 
Yij and dYij are defined). 

To define and dY^, we first select a set of points 
as follow: in each set i 1 ^ i, we select two points j 1 and 
j' + 1 such that Xi'ji < Xij < XiUj'+n, if there are no 
such points in a set, we do not select any point from 
that set (the set does not determine the position of the 
master curve for this value of x) . If this procedure selects 
no point at all then and dY^ are undefined for point 
ij and it does not contribute to S (this happens if set 
i is alone in this region of x and is the master curve 
by itself). We now compute the linear fit through the 
selected points (xi, yi, dyi), I = 1 . . .m and Yij is the value 
of that straight line at x^ and dY^ is the associated 
standard error, namely 



V- — 
1 y ~ 



(A2) 



and 



dY? = - (Kxx - 2 Xij K x + 4 K) (A3) 

with wi = l/dyf, K = ^2w h K x = Y, w i x u K y = 

Yj w iVu K xx = J2 w i x t> K xy = J2 w i x iVi and A = 
KK — K 2 
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The quality S measures the mean square distance to should thus be around one if the data really collapse to 
the master curve of the sets in unit of standard errors. It a single curve and much larger otherwise. 
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